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We theoretically study the atomic structure and energetics of silicon and silicon-nitrogen impuri¬ 
ties in graphene. Using density-functional theory, we get insight into the atomic structures of the 
impurities, evaluate their formation energies and assess their abundance in realistic samples. We 
find that nitrogen, as well as oxygen and hydrogen, are trapped at silicon impurities, considerably 
altering the electronic properties of the system. Furthermore, we show that nitrogen doping can 
induce local magnetic moments resulting in spin-dependent transport properties,even though nei¬ 
ther silicon nor nitrogen impurities are magnetic by themselves. To simulate large systems with 
many randomly distributed impurities, we derive tight-binding models that describe the effects 
of the impurities on graphene tt electron structure. Then by using the linear-scaling real-space 
Kubo-Greenwood method, we evaluate the transport properties of large-scale systems with random 
distribution of impurities, and find the fingerprint-like scattering cross sections for each impurity 
type. The transport properties vary widely, and our results indicate that some of the impurities can 
even induce strong localization in realistic graphene samples. 

PACS numbers: 61.48.Gh, 61.72.-y 


I. INTRODUCTION 

Graphene is a remarkable two-dimensional (2D) mate¬ 
rial due to its unique mechanical and electronic prop¬ 
erties, which makes it a good candidate for modern 
nanoscale devices and applications. While pristine 
graphene is a semimetal with a high carrier mobility^"— , 
for some applications it would be desirable to open a 
band gap in it, for example by cutting it to ribbons^”— or 
introducing impurities and defects^i'— . Moreover, impu¬ 
rities are often the dominant scatterers that control the 
intrinsic electronic and transport properties of realistic 
systems. 

Silicon is commonly present in nature, so that Si im¬ 
purities can appear in synthetic graphene during the 
growth processes or be introduced later when graphene 
is used together with standard silicon-based electronics 
components. It is even preferred to grow graphene on 
silicon wafers as one does not have to transfer it else¬ 
where after growth. Specifically, graphene grown at high 
temperatures by chemical vapor deposition can intro¬ 
duce silicon and oxygen impurities originating from the 
quartz (Si02) substrate or the apparatus itself. Also us¬ 
ing silicon carbide (SiC) ^^i^^ to grow graphene can pro¬ 
duce silicon impurites. Another avenue would be to de¬ 
liberately introduce silicon impurities by post-synthesis 
treatmentsi^ii^, such as by low-energy ion irradiation, 
similar to direct ion implantation of N and B atoms into 
single graphene sheetsi^, or deposition of Si atoms on 
ion or electron-beam treated graphene with irradiation- 
induced vacancies. Silicon atoms filling monovacancies 
and divacancies in graphene have been explicitly iden¬ 
tified in experimentsi^ii^. Their formation energies are 
expected to be fairly low compared to other period 3 el¬ 
ement substitution defectsi^, and they are stable enough 


for the electron beam not to break the atomic structure 
or the bonding easily^S. 

Si impurities can further pick-up various atoms from 
the environment. Recent experimental studies have re¬ 
ported the presence of individual defects formed by co¬ 
occurring silicon and nitrogen impuritie o^^i^^ . Zhou et 
al. showed that surface plasmons are locally enhanced at 
both silicon and a silicon-nitrogen impurities^S. There¬ 
fore, such silicon impurities occurring in graphene could 
in principle be used as plasmonic waveguides. Addition¬ 
ally, they could be useful in optoelectronic devices^l. Ni¬ 
trogen impurities have been extensively studied, as they 
dope their surroundings in graphen o^^i^^ , but the role 
of nitrogen binding to the silicon impurities, forming 
silicon-nitrogen defects, is not fully understood. More¬ 
over, the presence of oxygen can also affect the formation 
of silicon impurities. We answer these issues by finding 
the defect geometries and formation energies for vari¬ 
ous silicon, silicon-oxygen and silicon-nitrogen impurities 
through comprehensive density-functional-theory calcu¬ 
lations in a supercell geometry. 

Such impurities are also interesting in the context of 
magnetism. Local magnetic moments can be created in 
graphene in several ways. Examples include graphene 
nanoribbon zigzag edges^l, monovacancies and hydrogen 
adatoms2^, and transition metal substitutions^^. Local 
magnetic impurities interact strongly with the conduc¬ 
tion electrons of graphene, as has been shown by mea¬ 
suring the Kondo effecfj^i. We demonstrate that many of 
the nitrogen-doped silicon impurities exhibit finite spin 
moments. Such defects can be important in graphene- 
based spintronics applications, and we evaluate the spin- 
dependent electronic and transport properties for the 
most stable defect types. 

Electronic transport in systems with silicon point- 
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defects has been studied in a ribbon geometry by Lopez- 
Bezanilla et They found that the formation en¬ 

ergies of the silicon point-defects are smaller closer to 
the ribbon edges, and the transmission functions for 
graphene nanoribbons with silicon defects at the edges 
were computed. Furthermore, Cheng et al^^ studied the 
electronic and transport properties of the SiN^, defects 
in armchair nanoribbons. However, it is unclear what 
are the transport properties of realistic two-dimensional 
graphene systems with numerous randomly positioned 
defects. We simulate transport in such a realistic set¬ 
ting, and find the characteristic fingerprint-like scattering 
cross sections for each defect type. 

The paper is organized as follows. In Sections II A-D, 
we present first principles results for the geometries and 
formation energies of silicon, silicon-oxygen, and silicon- 
nitrogen impurities, and also touch upon the effects of 
hydrogen adsorption on them. In section II E, we eval¬ 
uate the electronic band structures and density of states 
of systems containing silicon and silicon-nitrogen impu¬ 
rities. In Section III A-B we derive tight-binding models 
to describe the effects of the impurities on the electronic 
structure. In Section III C, we compare the density of 
states of the system with impurities calculated within the 
periodic supercell approach with that of many randomly 
placed impurities in large realistic systems. In Section 
HID, we report the results of a comprehensive real-space 
Kubo-Greenwood study of the electronic transport of sili¬ 
con and silicon-nitrogen impurities, where we also discuss 
the localization effects in these systems. 


II. FIRST PRINCIPLES CALCULATIONS 

First principles calculations were performed by the all¬ 
electron density-functional theory package FHI-aims^S. 
The code uses local basis functions specified for each 
atom, and we have chosen the default tight basis sets pro¬ 
vided in the package. As the exchange-correlation energy 
functional we used the generalized-gradient approxima¬ 
tion as parametrized by Perdew, Burke, and Ernzerhof 
(PBE)^. 

The defects were placed in periodic supercells of 8 x 8 
and 11 X 11 graphene unit cells. Without any defects 
they would contain 2x8x8= 128 and 2x11x11 = 242 
carbon atoms, respectively. The lattice vectors of the su¬ 
percells were fixed to those of pristine graphene (with a 
lattice constant of 2.467 A in a two-atom unit cell), even 
after the point-defects had been introduced. The atom 
positions were relaxed until the forces between the atoms 
were smaller than 10“^ eV/A. The self-consistency cycle 
was considered converged if the change in the volume- 
integrated root-mean square of the charge density, and 
changes in the sum of eigenvalues and total energy were 
below 10“^ [electrons], 10“^ eV, and 10“^ eV, respec¬ 
tively. A Monkhorst-Pack grid of 8 x 8 x 1 k-points was 
employed in the total energy calculations, and the band 
structure and density of states (DOS) calculations uti¬ 


lized a F-centered grid with at least 15 x 15 x 1 k-points 
in order to achieve a better coverage of the first Brillouin 
zone. 

Performing the calculations with two distinct supercell 
sizes lets us explicitly test whether the results have con¬ 
verged as a function of system size. Systematic conver¬ 
gence studies have been done for the silicon substitution 
defect in graphene in Ref. [T2|, where convergent results 
were obtained already with a rather small 4x4 super¬ 
cell, as contrary to vacancies^, substitutional impuri¬ 
ties and adatoms do not give rise to long-ranged strain 
fields. Studying structural point-defects in graphene, a 
7x7 supercell is typically large enough for geometry 
relaxation^. However, to obtain convergent electronic 
properties and spectra, a larger supercell is typically 
needed. Specifically, for isolated nitrogen substitution 
defect, or any other defect with large effective on-site 
potential at the impurity site, the supercell needs to be 
much larger—. The same is true for the adequate de¬ 
scription of magnetic properties of some defects24. 

To estimate the likelihood of a defect to form, and to 
compare the energies of different defects in graphene, we 
define the formation energy of a defect as 

-£/f E ^ ^ (l) 

where E is the total energy of a supercell with the defect 
in question, i sums over atom types, rii is the number 
of atoms of type i, and /i^ is the chemical potential of 
atom type i. Now, the carbon chemical potential //c was 
chosen as the energy of graphene per carbon atom, so 
that defect-free graphene has zero formation energy. For 
the silicon chemical potential //si, we chose the energy 
per atom of crystalline silicon in a diamond cubic lattice. 
The hydrogen, nitrogen, and oxygen chemical potentials 
/iH, /^N, and /io were chosen as half the energies of the 
H 2 , N 2 , and O 2 molecules, respectively. The formation 
energies can be related to binding energies given the for¬ 
mation energies of isolated H, C, N, O, and Si atoms, 
that are 2.27 eV, 8.01 eV, 5.28 eV, 2.92 eV, and 4.66 eV, 
respectively. 

We would like to stress here that even though concen¬ 
tration n of a particular type of defects at temperature 
T can be evaluated as n ^ exp(—/3F^f) for the system 
in the thermodynamic equilibrium only, Eq. m makes it 
possible to qualitatively assess relative concentrations of 
different types of defects in the system subjected to mild 
chemical treatment, and even though in some cases after 
harsh treatment like irradiation, provided that the sys¬ 
tem was annealed after that. All defects were assumed 
to be neutral, as graphene is a semimetal. 

The studied defect cases and shorthand notations for 
labelling them are presented in Eig.[Tl When a single Si, 
N or O atom fills a monovacancy (IVac) in the graphene 
lattice, the substitutions are denoted as ISi, IN and 10, 
and when it fills a divacancy (2Vac), the substitutions 
are denoted as 2Si, 2N and 20, respectively. A composite 
defect consisting of two defects, for instance two substi- 
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FIG. 1. (Color online). Atomic structures, top and side views, for defects in graphene containing silicon (large green spheres), 
oxygen (red spheres), nitrogen (blue spheres), and hydrogen (small yellow spheres) atoms. The geometries are relaxed in the 
11 X 11 supercells, and only the defect neighbourhood is shown, (a)-(e) Silicon impurities, (f)-(j) Silicon-oxygen impurities 
(k)-(o) Silicon-nitrogen impurities, (p)-(t) Hydrogen adatoms on silicon and silicon-nitrogen impurities. 


tutions, located at neighboring sites is denoted as lSi-10. For adatoms on pure graphene, we specify the absorption 
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site like Si-top, and further when the adatom is on top of 
a defect, we combine the notation to form 0-top@lSi for 
a oxygen adatom on top of a silicon in a monovacancy. 
Below, we first show the results for the defects containing 
only silicon and oxygen, see Figs. [Il(a)-(j). After that, 
nitrogen is added, see Figs. [TJ(k)-(o). Finally, the hydro¬ 
gen adatom cases shown in Figs. [Tl(p)-(t) are studied. 


A. Silicon and oxygen impurities 

Silicon and oxygen impurities can appear in graphene 
for example due to contamination in growth processes 
where Si02 or SiC are present, or due to purposeful 
post-synthesis treatment. To determine which impuri¬ 
ties are the most likely to occur, we relaxed various silicon 
and oxygen point-defect atomic geometries and evaluated 
their formation energies. 

We first analyze the results for the silicon defects. The 
ISi defect, where a single carbon atom is substituted by a 
heavier silicon atom, prefers the tetrahedral sp^ type hy¬ 
bridization, resulting in an out-of-plane silicon position, 
as seen in Fig. [TJa). Consequently, the surrounding car¬ 
bon atoms move accordingly, and the resulting curvature 
takes several bond lengths to relax, making it effectively 
a slightly extended defect. On the other hand, the silicon 
atom in a 2Si defect, as shown in Fig.lTJb), forms sp^d hy¬ 
brid orbitals that bond with the four neighboring carbon 
atoms, and the geometry is only slightly perturbed in the 
out-of-plane direction. The defect geometries and orbital 
hybridizations of the ISi and 2Si defects have been estab¬ 
lished both by c omp utational and experimental methods 
in Refs, pisl andfl^. 

The defect formation energies are listed in Table [H 
The ISi and 2Si defect formation energies are 3.77 eV 
and 4.57 eV, respectively, based on the 11 x 11 super¬ 
cell calculation. Increasing the supercell size typically 
lowers the formation energy, except in some cases where 
the computational supercell is too small, such that there 
is not enough space for bulk graphene between the de¬ 
fects, or due to other finite size effects^. The chosen 
chemical potentials affect the absolute values of the 
formation energies, but when comparing the formation 
energies of defects with the same number of atoms 
the riiiiii terms cancel in Eq. (HD. Therefore, regard¬ 
less whether the silicon atom migrates from the sub¬ 
strate or arrives as part of a molecule exposed to the 
surface, as long as the source of the silicon remains the 
same, 2Si has a higher formation energy than ISi by 
F;f[2Si]-F;f[lSi] = E[2Si]-E[lSi]^fic ^ 0.8 eV. The left¬ 
over carbon atoms are assumed to form bulk graphene, 
but if they obtained an energy higher than jic^ ISi would 
become even more favorable compared to 2Si. On the 
other hand, if the defects were to form on the spot, the 
energy per ejected carbon atom is lower for 2Si. In any 
case, the sp^ hybridization in ISi with three-fold coordi¬ 
nation is energetically preferred to the sp^d hybridization 
in 2Si with four-fold coordination, even if the ISi defect 


causes a local curvature of the graphene sheet. 

A divacancy can be relaxed further by a bond rotation 
to obtain a 555-777 defect with even lower energy^. We 
studied the case where a silicon atom substitutes one of 
the carbon atoms at the 555-777 defect. It turned out 
that silicon atom prefers to substitute not the middle site 
but a site neighboring it, as shown in Fig.^c). However, 
the formation energy of such a defect is much higher than 
that of ISi and 2Si. 

A silicon adatom on pristine graphene takes a posi¬ 
tion on top of a carbon atom (Si-top), but it is slightly 
shifted towards the hexagon of the graphene backbone, 
see Fig. [TJd). Its formation energy is 4.33 eV, being be¬ 
tween the ISi and 2Si formation energies. It is only barely 
lower than the formation energy of an isolated silicon 
atom that is 4.66 eV. Moreover, the total energy of the 
silicon adatom in a bridge position on top of a carbon- 
carbon bond (Si-bridge), shown in Fig.lTJe), is only a few 
meV higher, further indicating that silicon adatoms are 
rather mobile on graphene. This, in turn, implies that 
silicon adatoms are easily trapped by monovacancies and 
divacancies, because the formation energies of ISi or 2Si 
are much lower than the sum of the formation energies 
of an isolated silicon atom or Si-top/bridge, and I Vac 
or 2Vac. The energy gain is roughly 7-8 eV. Further¬ 
more, a system consisting of a pair of neighbouring ISi 
defects (I Si-1 Si), or a pair of neighbouring ISi and 2 Si 
defects (ISi-2Si), has a total energy that is 1.6 eV, or 
2.3 eV, lower than the energy of two fully separated de¬ 
fects, respectively. Such composite defects minimize the 
total distortion, such as the curvature of the graphene 
lattice, and therefore reduce the total energy. However, 
silicon substitution defects are not expected to be mobile 
(at least under ambient conditions) after they have been 
formed. 

Oxygen could also play an important role in forming 
impurities in graphene. It is ubiquitous, and it is com¬ 
monly present during the graphene growth process or 
measurements especially on a Si02 substrate. Oxygen in 
graphene prefers the bridged adatom position (0-bridge, 
not shown) as opposed to filling a monovacancy or a di¬ 
vacancy, 10 or 20 (not shown). Namely, the formation 
energy of an 0-bridge is about 2 —3 eV smaller than those 
of 10 and 20, and also about 0.8 eV smaller than the 
case where the oxygen adatom is located exactly on top 
of a carbon atom. Contrary to silicon, an oxygen atom is 
too small to completely fill a divacancy, and a 20 defect is 
better described as 10-1 Vac. Curiously, two neighbour¬ 
ing oxygen substitutions (10-10), shown in Fig.lTJj), has 
a remarkably low formation energy of — 1.12 eV. However, 
there can be a high energy barrier to reach such a con¬ 
figuration, as hinted by the high formation energy of 20 
(10-1 Vac) that is 2.96 eV. 

It is not a surprise that, when considering compos¬ 
ite defects containing both silicon and oxygen, silicon 
prefers to form ISi and 2Si defects, and oxygen prefers 
the adatom position. However, the lowest formation en¬ 
ergy of 1.40 eV is obtained for an oxygen adatom located 
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TABLE L Defect spin moments and formation energies. The spin moments are for relaxed 11x11 supercells, and the formation 
energies are also given for the 8x8 supercell. Energies are in units of eV. The spin moment values inside brackets are 
atom-projected spin moments using the Hirshfeld analysis. 


Impurity atoms 

Defect 

Spin moment { i ^ b ) 

£5(8 X 8) £5(11 

X 11) £;f(Ref.) 

- 

IVac 

1.44 [0.69(C)] 

7.68 

7.65 

7.5S 


2Vac 

- 

7.56 

7.36 

335 

Si 

Isolated Si 

2.0 

4.66 

4.66 



ISi 

- 

3.79 

3.77 



2Si 

- 

4.60 

4.57 

4.41^ 


ISi @ 555-777 

- 

6.89 

6.96 



Si-top/bridge 

- 

4.33 

4.33 



ISi-lSi 

- 

5.99 

5.97 



lSi-2Si 

- 

6.09 

6.06 


0 

Isolated 0 

2.0 

2.92 

2.92 



10 

- 

2.43 

2.44 



20 (lO-lVac) 

- 

3.03 

2.96 



0-bridge 

- 

0.43 

0.43 



10-10 

- 

-1.10 

-1.12 


Si + 0 

lSi-10 

- 

2.78 

2.77 



2Si-10 

- 

2.26 

2.22 



lSi-20 

- 

5.03 

4.98 



2Si-20 

- 

4.66 

4.52 



0-top @ ISi 

- 

1.43 

1.40 



0-bridge @ 2Si 

- 

1.69 

1.66 



Si-bridge @10 


4.06 

4.04 


N 

Isolated N 

3.0 

5.28 

5.28 



IN 

- 

0.93 

0.91 

0.92a 


2N (IN-lVac) 

- 

5.52 

5.45 



N-bridge 

0.72 [0.53(N)] 

4.38 

4.37 



IN • • • IN (3rd-NN) 

- 

1.91 

1.91 


Si + N 

ISi-lN 

1.00 [0.44(Si), 0.07(N)] 

3.34 

3.33 



2Si-lN 

- 

3.25 

3.22 

2.59^ 


lSi-2N 

- 

8.26 

8.21 



2Si-2N 


7.15 

7.14 



N-bridge @ ISi 

- 

5.79 

5.75 



N-bridge @ 2Si 

0.67 [0.02(Si), 0.26(N)] 

5.87 

5.83 



Si-bridge @ IN 

1.00 [0.74(Si), 0.02(N)] 

5.16 

5.16 


H 

Isolated H 

1.0 

2.27 

2.27 



H-top 

1.00 [0.04(H)] 

1.45 

1.45 


Si + H 

H-top @ ISi 

- 

3.23 

3.21 



H-top @ 2Si 

- 

4.61 

4.56 


N + H 

H-top @ IN 

- 

2.68 

2.68 


Si + N + H 

H-top(Si) @ ISi-lN 

- 

2.03 

2.02 



H-top(N) @ ISi-lN 

- 

4.54 

4.52 



H-top(Si) @ 2Si-lN 

- 

2.91 

2.88 



H-top(N) @ 2Si-lN 

- 

4.12 

4.08 



on top of a ISi defect (0-top @ ISi), shown in Fig. [T]^h), 
and not in a bridged position. The oxygen adatom bonds 
with the silicon atom in ISi much more preferably than 
with graphene, and surprisingly, also more preferably 


than forming an O 2 molecule. An oxygen bridge located 
on a 2Si defect (0-bridge @ 2Si), shown in Fig.^i), also 
has an unexpectedly low formation energy of 1.69 eV. 
These results indicate that oxygen is captured at the sil- 
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(a) 1Si-1N 
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FIG. 2. (Color online). The total energies E and spin mo¬ 
ments for the (a) ISi-lN and (b) 2Si-lN defects in relaxed 
11 X 11 supercells, where the N atom is separated to lattice 
sites further away from the Si atom. The zero energy val¬ 
ues correspond to the nearest-neighbour cases of the ISi-lN 
defect in (a) and the 2Si-lN defect in (b). The Si and N 
atom-projected spin moments are evaluated using the Hirsh- 
feld analysis. 


icon impurities. 

When embedding oxygen in the graphene lattice, the 
composite defects such as lSi-10 and 2Si-10, shown in 
Fig. [1] (f) and (g), have comparably small formation en¬ 
ergies as well, 2.78 eV and 2.26 eV, respectively. Thus it 
costs less energy to form oxygen substitutions next to ISi 
or 2Si than in pristine graphene. There are many possible 
configurations that could be called 1 Si-20 or 2Si-20, but 
they all have systematically at least about 2 eV higher 
energies than lSi-10 and 2Si-10. The same goes for re¬ 
versing the roles of silicon and oxygen, namely a defect 
consisting of a silicon adatom located near a 10 defect 
has a much larger formation energy. Otherwise, the sili¬ 
con and oxygen pair well to form composite point-defects. 


B. Nitrogen doping and finite spin moments 

A nitrogen impurity in graphene dopes its surround¬ 
ings, which can be used to tune the electronic properties 
of graphene. A nitrogen substitution defect (IN), shown 
in Fig. [D^o), is energetically by far the most stable de¬ 
fect configuration that contains nitrogen only. The differ¬ 
ences in formation energy to a bridged nitrogen adatom 
(N-bridge, not shown) and to a divacancy substitution 
(2N or IN-lVac, not shown) are 3.4 eV and 4.6 eV, re¬ 
spectively. Since a nitrogen atom has one more electron 
than a carbon atom, and a substitutional nitrogen atom 
forms sp^ hybridized bonds in graphene, some of the elec¬ 
tron density is donated to the nearby atoms, effectively 
doping the neighbour hood^^. Thus, we are mostly inter¬ 
ested in the IN-doped silicon impurities, to find how the 


added nitrogen impurity changes the properties of the 
silicon point-defects. 

The straightforward case, with neighboring silicon and 
nitrogen substitutions (ISi-lN), shown in Fig. [TJk), has 
a comparably low formation energy of 3.33 eV. Inter¬ 
estingly, even if neither ISi or IN is magnetic individ¬ 
ually, the ISi-IN defect has a finite total spin moment 
of l.OO/iB* Table [I] also shows the atom-projected Hir- 
shfeld spin moments, indicating that almost half of the 
total spin density resides near the silicon atom. This is 
explained by that the odd electron originating from the 
nitrogen atom cannot pair and forms a spin-polarized 
impurity state centered around the silicon atom. The 
nitrogen doping does not distort the atomic structure es¬ 
pecially much compared to the ISi case, suggesting that 
the silicon still hybridizes as sp^. In fact, as the nitro¬ 
gen atom is able to form sp^ bonds, the ISi-IN defect is 
comparable to the ISi defect with the distinction of the 
nitrogen doping. Thus, in a sense, the IN defect dopes 
the ISi defect and creates a magnetic moment. However, 
the difference in the total energies of the spin-polarized 
and spin-unpolarized solutions is rather small, 72 meV 
in the 11 x 11 supercell, meaning that one has to reach 
low temperatures not to mix these states. 

The nitrogen doped 2Si defect, namely 2Si-lN shown in 
Fig. mi), does not have a finite spin moment. However, 
its formation energy is extremely low, 3.22 eV, which 
is 0.11 eV lower than that of the ISi-IN defect. Fur¬ 
thermore, the 2Si-lN defect is relaxed perfectly in-plane, 
unlike the 2Si defect. Both the ISi-IN and 2Si-IN de¬ 
fects have been observed in experiment o^^i^^ , reflecting 
the fact that their formation energies are low in compar¬ 
ison to any competing defects. In addition, Zhou et al3 
evaluated the energies required to remove either the sili¬ 
con or the nitrogen atom from the 2Si-IN configuration, 
obtaining the high values of 7.45 eV and 7.01 eV, re¬ 
spectively, which further highlights the stability of these 
defects. 

Besides the ISi-IN defect, some of the nitrogen- 
doped defects are magnetic, such as the bridged nitro¬ 
gen adatoms on graphene (N-bridge, not shown) and 
on 2Si (N-bridge @ 2Si, not shown), and even a silicon 
bridge on IN (Si-bridge @ IN), shown in Fig. Cu¬ 

riously, a nitrogen bridge on ISi (N-bridge @ ISi), shown 
in Fig. m^m), is not magnetic even if a nitrogen bridge 
on graphene is. However, all their formation energies are 
higher than those of ISi-IN and 2Si-IN, but still lower 
than having the adatoms moved on graphene. In this 
sense, nitrogen and silicon adatoms are attracted to ISi, 
2Si and IN impurities, and then they are possibly relaxed 
to ISi-IN and 2Si-IN defects. 

We further tested both the ISi-IN and 2Si-lN defects 
by separating the silicon and nitrogen atoms from the 
nearest-neighbor sites. In practice, the nitrogen atom 
was swapped with a carbon atom farther away in the 
lattice, after which the geometry was fully relaxed. The 
resulting total energies are shown in Fig. [2] (top) as a 
function of the distance between the silicon and nitrogen 












7 


atoms. For both ISi-lN and 2Si-lN, separating the sili¬ 
con and nitrogen atoms costs more than 1 eV in energy. 
Intuitively, the ISi and IN defects being in the nearest- 
neighbor sites causes the least amount of distortion in 
the honeycomb lattice and there is only minimal associ¬ 
ated energy cost. This also validates the nearest-neighbor 
configurations ISi-lN and 2Si-lN as the most stable, and 
therefore the most important to study. In the large Si-N 
distance limit, the energies in Fig. [2] approach the for¬ 
mation energies of separated ISi and IN, and 2Si and 
IN defects, namely the values E'f[lSi-lN] -1-1.35 eV and 
F^f[2Si-lN] -h 2.26 eV, respectively. In this sense, ISi and 
2Si defects separated by more than about 3A from a IN 
defect can already be considered almost fully separated 
in terms of the energetics. 

The total and atom-projected Hirshfeld spin moments 
for the defect configurations with the separated silicon 
and nitrogen defects are shown in Fig. [2] (bottom). The 
systems with ISi and IN defects are magnetic practi¬ 
cally only if the silicon and nitrogen atoms reside as the 
nearest- or third-nearest neighbors. On the other hand, 
the 2Si and IN defect configurations obtain finite spin 
moments only if the silicon and nitrogen atoms are far¬ 
ther than about 4A away from each other. However, as 
the separation distance grows, it is expected that the 
system becomes spin-unpolarized, as in the case of indi¬ 
vidual 2Si and IN defects. 

Since the ISi-IN defects are magnetic, it is meaningful 
to ask how the defect spin moments interact and align in 
a system with many such defects. To study this, we put 
two ISi-IN defects in a 11 x 11 super cell such that the 
defects formed a honeycomb superlattice of their own. 
Based on calculations using the default light basis for 
each atom, due to the large number of possible configu¬ 
rations, it turned out that the ISi-IN defects prefer an¬ 
tiferromagnetic ordering, regardless of which sublattices 
the impurity atoms reside on, and whether the defect out- 
of-plane positions are in the same direction or not. The 
antiferromagnetic ordering is consistently tens of meV 
lower in energy compared to the ferromagnetic or spin- 
unpolarized solutions. The energy differences are rather 
small, and the spin moments can be assumed to interact 
only weakly especially in low defect concentrations. 


C. Hydrogen adatoms on silicon impurities 

Interestingly, hydrogen adatoms located on top of sili¬ 
con and silicon-nitrogen impurities have remarkably low 
formation energies, see Table [H Hydrogen adatoms on 
graphene (H-top), shown in Fig. mt), have an energy 
penalty of 1.45 eV per hydrogen atom, when using a H 2 
molecule as a reference for the chemical potential fin in 
Eq. ([1]). However, the formation energies of silicon impu¬ 
rities actually decrease when a hydrogen adatom is added 
on top the silicon impurity. The decrease in energy is the 
most notable for the ISi and ISi-IN defects, being 0.56 
eV and 1.31 eV, respectively. The IN-doping plays a sig- 



FIG. 3. The first Brillouin zone of the computational super¬ 
cells. The reciprocal primitive vectors are denoted as hi and 
&2. When evaluating the band structure, we focus on the tri¬ 
angular paths between nearby high symmetry points, namely 
FMAF, TM'K'T and TM"K"T. 

nificant role by lowering the energy even further, even 
though the hydrogen adatom is clearly located on top of 
the silicon atom, see Fig. Up) and (s). The low forma¬ 
tion energies overall hint that the silicon-nitrogen defects 
could be reactive not only to hydrogen but also to vari¬ 
ous molecules as well. This is the prerequisite for using 
the defects as sensors. 

A hydrogen adatom on graphene has a total spin- 
moment of 1.0//B- This follows from the fact that the 
hydrogen atom bonds to a carbon atom, passivating the 
carbon Pz orbital and effectively resulting in an empty 
site as seen from the point of view of the tt electron 
bands. By Lieb’s theoren>^, a finite spin moment will 
result from the sublattice imbalance. Quite surprisingly, 
magnetization does not occur when the hydrogen adatom 
is located on top of a ISi defect, or a ISi-IN defect. In 
fact, all the silicon and silicon-nitrogen impurities with a 
hydrogen adatom considered here, see Table [H are non¬ 
magnetic. Therefore hydrogen adatoms, which readily 
trap at the silicon and silicon-nitrogen impurities, passi¬ 
vate any defect-induced magnetism. 


D. Electronic band structures 

In the following, we studied the DFT band structures 
and DOS of systems with 8x8 unit cells of graphene that 
contain a single ISi, 2Si, ISi-lN, or 2Si-lN defect. This 
corresponds to an infinite, periodic array of the impuri¬ 
ties. The band energies are evaluated in the reciprocal 
space along the line segments between the nearby high 
symmetry points in the first Brillouin zone (BZ), as il¬ 
lustrated in Fig. [3l From the outset, we restricted the 
study only to one of the halves of the first BZ, since time 
reversal symmetry is not broken, and therefore inversion 
around the F point, namely k — k, does not alter state 
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FIG. 4. (Color online). Band structures, DOS and atom-projected pDOS of 8 x 8 supercells with (a) a ISi defect, and with 
a spin-polarized ISi-lN defect showing the (b) majority and (c) minority spin components separately. The DFT results are 
the solid lines, and the fitted tight-binding model results are the dotted lines. The pristine graphene band structure and DOS 
(thick grey lines) are shown in (a). The pDOS curves are the total atom-projected pDOS (black), the s-type pDOS (yellow), the 
p-type pDOS (light blue), and the d-type pDOS (red). The DOS and pDOS have been broadened using Gaussian broadening 
of 30 meV. Zero energy is at the Fermi energy, the charge neutrality point. The atom labels match those in Figs. [TJa) and (k). 
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(^) 2Si Qi m r.o r.Q r.A 



FIG. 5. (Color online). Band structures, DOS and atom-projected pDOS for (a) the 2Si and (b) the 2Si-lN defects. The 
DFT results are the solid lines, and the fitted tight-binding model results are the dotted lines. The pDOS plots show the total 
atom-projected pDOS (black), the s-type pDOS (yellow), the p-type pDOS (light blue), and the d-type pDOS (red). The 8x8 
supercell has been used, and DOS and pDOS have Gaussian broadening of 30 meV. Zero energy is the Fermi energy, the charge 
neutrality point. The atom labels match those in Fig. [TJb) and (1). 


energies. This also holds for the majority and minority 
spin channels in the spin-polarized case. 

The band structure, DOS and atom-projected partial 
DOS (pDOS) of the system with a ISi defect are shown 
in Fig. IHa). The band structure closely matches the 
pristine graphene band structure (thick grey line). In 
fact, one sees in the presented energy range the pristine 
graphene occupied and unoccupied tt electron bands that 
have been folded multiple times in order to match the 
8x8 supercell as a periodic unit. A supercell with a 
defect cannot be unfolded back into smaller periodic unit 
cells, as there are energy gaps between the bands at the 
high symmetry points. Curiously, with supercells of sizes 


p X p where p is divisible by three, the K and K' points 
are actually mapped onto the F point, and there would 
be no energy gap at However, in our case with 

p = 8, there is a small energy gap at the K point in the 
ISi band structure, but otherwise the Dirac cone seems 
to be unperturbed. 

Consequently, the ISi DOS is similar to the graphene 
DOS, indicated by the thick grey line in Fig.lU^a), except 
that there are additional states at energies around 0.8 
eV. The opening of energy gaps in the band structure 
results in double peaks in the DOS if the peak broad¬ 
ening is sufficiently small. But most importantly, the 
peaks in the DOS are in general explained by the flatter 
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energy bands, for instance resulting from localized impu¬ 
rity states that are only weakly k-dependent. However, 
substituting a carbon atom by a valence isoelectronic sil¬ 
icon atom does not introduce a new impurity band, but 
rather the tt electron bands become perturbed. Peaks 
are then formed as there are stronger resonances at cer¬ 
tain energies that are also expected to depend on the 
dimensions of the periodic array of defects. The pDOS 
around the ISi defect, in Fig. SJ^a), shows that the peaks 
in the density at around 0.8 eV are mostly localized at 
the silicon atom and the nearest and third-nearest car¬ 
bon atoms Cl and C3, see Fig. [TJa) for the atom labels. 
We can also see that carbon atoms exhibit only p-type 
orbital occupations as expected, but the silicon atom has 
some contribution from s- and even d-type orbitals. 

The ISi-lN band structure and DOS are shown in Fig. 
[T{b) and (c) for the majority and minority spin compo¬ 
nents from a calculation with collinear spin. The paths 
in the reciprocal space along which the band energies are 
shown now include the three triangular paths FMi^F, 
TM'K'r^ and rM"K"r due to breaking of the graphene 
lattice symmetries. There is a spin moment of 0.9SfiB 
that corresponds to an almost fully occupied band in the 
majority spin channel that is unoccupied in the minor¬ 
ity spin channel. The origin of the spin moment is the 
nitrogen atom that introduces the odd electron that is 
attracted to the silicon atom. As a consequence, there 
is an almost flat impurity band just below the Fermi en¬ 
ergy, and a clear bump in the DOS and atom-projected 
pDOS. The density is mostly localized at the silicon atom 
with some weight also at nearby atoms of the adjacent 
sublattice. In addition to the p-type orbital occupations, 
the impurity state also exhibits some s- and d-type or¬ 
bital occupation at the silicon atom, and small d-type 
occupations even at the second-nearest neighbor carbon 
atom C2. 

For both spin components, the band energies along the 
three paths TMKT, TM'KT and rM"K'T have only 
minor differences, suggesting that there is no strong di¬ 
rectional dependence. Therefore, the electronic states at 
the nitrogen atom are expected to be similar to the states 
at the carbon atoms that are nearest-neighbours of the 
silicon atom. In fact, the band structure of the minor¬ 
ity spin component, shown in Fig.lTJc), closely resembles 
the band structure of the ISi defect. The Dirac cone at 
the K-point is intact, being only slightly shifted to the 
occupied side due to the nitrogen doping. However, with 
the larger system of 11 x 11 graphene unit cells, the spin 
moment is already 1.00//^, and therefore the Dirac cone 
in the minority channel is expected at the Fermi energy. 
We can call such a material half-semimetal, since the ma¬ 
jority spin channel has only an impurity band just below 
the Fermi energy, and the band structure of the minor¬ 
ity spin channel is that of a semimetal. Furthermore, the 
main qualitative difference between the pDOS of periodic 
arrays of ISi and ISi-IN defects is that the silicon pDOS 
shows two peaks in the ISi-IN minority channel, whereas 
the silicon pDOS at the ISi defect has three peaks. More¬ 


over, there is no clear occupation of d-type orbitals in the 
minority channel of the ISi-IN defect system. 

The band structure, DOS and pDOS of a system with 
a 2Si defect are shown in Fig. [5fa). The breaking of the 
sublattice and the three-fold rotational symmetries have 
a profound effect on the band structure. Namely, there 
is a strong directional dependence on the wave vectors 
k, and even the Dirac cones are slightly shifted from the 
K, K' and K" points. Furthermore, there is a clear 
bump in the DOS roughly at 0.2 eV. This shows in the 
pDOS as localized density from p- and d-type orbitals at 
the silicon atom and as density from p-type orbitals at 
the nearby carbon atoms. Interestingly, there are states 
in the occupied side that do not have much density at 
the silicon atom, but they do have density at the nearby 
carbon atoms. Therefore there are also impurity states 
surrounding the silicon site. 

The band structure of a system with a 2Si-lN defect, 
shown ininfb), is necessarily doped due to the odd elec¬ 
tron from the nitrogen atom. There is a band at the 
Fermi energy that is half filled, and as a result the conical 
dispersions are shifted roughly by 0.4 eV to the unoccu¬ 
pied side. In fact, the 2Si-lN band structure looks more 
like doped ISi band structure than 2Si band structure. 
Compared to 2Si band structure, where the bands have 
strong directional dependency, the valence bands in the 
2Si-IN band structure are almost flat resulting in many 
distinct peaks in the DOS and pDOS in the defect neigh¬ 
borhood. These states are again not only localized im¬ 
purity states at the silicon but also at the nearby carbon 
atoms, especially the nearest-neighbour carbon atom Cl. 
Curiously, there is an almost perfectly flat band at —1.26 
eV, implying a state localized at the defect, and that the 
related orbitals have a symmetry other than that of the 
TT electrons. 


III. TIGHT-BINDING RESULTS 


A. Tight-binding models 


We constructed tight-binding models for the silicon 
and silicon-nitrogen point-defects on graphene to further 
study their electronic and transport properties. This way 
we could simulate larger systems with random defect con¬ 
figurations, which is beyond typical first principles ap¬ 
proaches. 

The tight-binding model with an orthogonal state per 
each atom site, and simple hoppings between the sites, 
describes the graphene tt and tt* electron bands remark¬ 
ably well^i^. The tight-binding Hamiltonian is written 
as 


H = c\ci + 

i 


tij cjcj -h H.c. 

- ^<J 


( 2 ) 


where cj (q) creates (annihilates) a fermion in a state 
i, and H.c. stands for Hermitian conjugate. Here, the 
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coefficients Si describe on-site potentials, and tij are hop¬ 
pings between sites. We denote as the nearest-, 

second-nearest-, and third-nearest neighbor hoppings, re¬ 
spectively. 

The tight-binding models for systems containing de¬ 
fects were found by modifying the hoppings and on-site 
potentials in the vicinity of the defects, and fitting the 
outcome to the first principles band structures. Further¬ 
more, we checked that the resulting tight-binding local 
density of states (LDOS) matched the first principles 
atom-projected partial DOS around the Fermi energy, 
at least qualitatively. We used the 8x8 supercells that 
have a single defect each. The band structure fitting was 
performed in the first Brillouin zone in a uniform grid of 
50 X 50 X 1 k-points. We took 12 bands into the fitting 
procedure, the six highest valence bands and the six low¬ 
est conduction bands, determined at the F point, except 
for the 2Si-IN defect, where only four valence and four 
conduction bands were taken in order to avoid the flat 
band at —1.26 eV. 


B. Tight-binding parametrization 

First the carbon-carbon (C-C) hoppings ti,t 2 ,t 3 were 
obtained by fitting to the pristine graphene band struc¬ 
ture. We allowed ti and t 2 to be free fitting param¬ 
eters to scale the unit of energy and to freely break 
the electron-hole symmetry by tuning ^ 2 - We included 
the third-nearest neighbor hoppings ts by assuming that 
ts = (0.18/2.7), where the factor is motivated by previ¬ 

ous graphene tight-binding models^. If ts was not con¬ 
strained, altering ti and ts simultaneously would produce 
a continuous set of models with good fits. The resulting 
C-C hoppings are = (—2.855,—0.185,—0.190) 

eV, which were used for all carbon-carbon hoppings in 
our tight-binding calculations. 

We expect that the properties of both the silicon and 
silicon-nitrogen defects can be reproduced around the 
Fermi energy using simple tight-binding models. Namely, 
silicon and carbon atoms are valence isoelectronic, and 
nitrogen is also able to form sp^ bonds. The relevant tt 
and TT* electron bands are then perturbed by the defects 
in the effective tight-binding description. Moreover, we 
chose a minimal set of tight-binding parameters to en¬ 
hance the physical significance of each parameter. We 
took only the nearest-neighbor hoppings between the Si 
and C atoms, and between the Si and N atoms, to be 
nonzero. The hoppings between C and N atoms were 
taken identical to those of the C-C hoppings, which is 
similar to the nitrogen substitution model in Ref. [2^. 
The C on-site potentials were set to zero, whereas the Si 
and N on-site potentials could take non-zero values. 

With the parameter constraints and fitting method de¬ 
scribed above, the resulting tight-binding parameters for 
the ISi, ISi-lN, 2Si, and 2Si-lN defects are listed in Ta¬ 
ble El We provide spin-unpolarized and spin-polarized 
parametrization of the ISi-lN defect with distinct models 


TABLE II. Tight-binding parameters for the ISi, ISi-lN, 2Si 
and 2Si-lN defects. The C-N hoppings are taken as the same 
as the C-C hoppings, and the C on-site potentials are assumed 
as zero. The values are in units of eV. 


Defect Spin polarization 

Hopping ti\ 

Potential e: 



C-Si Si-N 

Si 

N 

ISi 

no 

-1.123 

0.118 


ISi-lN 

majority 

-0.776 -0.919 

-0.539 

-4.874 


minority 

-0.921 -1.483 

0.165 

-5.182 


no 

-0.819 -1.419 

-0.351 

-5.393 

2 Si 

no 

-0.849 

2.087 


2Si-lN 

no 

-1.627 -0.894 

2.722 

-2.833 


for the majority and minority spin components. Compar¬ 
ing to first principles, the tight-binding models produced 
remarkably accurate band structures, DOS, and LDOS, 
as shown by the dotted lines in Figs. 2] and El The 
tight-binding Fermi energies are at 0.537 eV, 0.593 eV, 
0.392 eV, and 0.274 eV for the systems with a ISi, ISi- 
IN (spin-polarized), 2Si, and 2 Si-lN defect, respectively. 
Below, we use these values as E" = 0 point, even if the 
defect concentration changes. 

The most visible signature of the ISi defect is the open¬ 
ing of the energy gaps between the lowest and second low¬ 
est conduction bands at the F point, as well as between 
the second and third lowest conduction bands at the M 
and K points. This can be explained by the weaker Si-C 
hopping compared to the C-C hopping value. An on-site 
potential at the silicon atom also breaks the translational 
symmetry, resulting in small energy gaps at many band 
crossings. In the case of the 2 Si defect, the C-Si hopping 
is even weaker, and we set a repulsive on-site potential 
at the silicon atom, since the sp^d hybridization takes 
more electrons to form bonds with the four neighboring 
carbon atoms. For the defects containing nitrogen, we 
used a rather low on-site potential at the nitrogen sites, 
because a nitrogen substitution dopes its surroundings 
by sharing some of its electron density. For the 2Si-IN 
defect, we did not include the state corresponding to the 
extremely flat band in the model. 

Curiously, a ISi-lN defect is spin-polarized. We ob¬ 
tained the tight-binding Hamiltonians for both the ma¬ 
jority, El, and the minority, E 25 spin components indi¬ 
vidually. To explicitly include the electron spin in the 
Hamiltonian, we can write 

H = Ei(g)|cri)(cri|+E 2 ( 8 )|cr 2 )(cr 2 | = E( 8 ) 1 + , (3) 

where the electron spin is first projected onto the defect 
majority and minority spin states, \(Ti) and |cr 2 ), and we 
have defined E = (Ei + E 2 )/ 2 , A = (Ei - E 2 )/ 2 , and 
= |cri)(<^i| — \(^ 2 ){cF 2 \‘ Due to the last term in the 
Hamiltonian, the electron spin can flip, if it is not aligned 
with the defect spin moment. The strength of the spin 
flipping is proportional to the difference of the majority 
and minority tight-binding parameters of A. In our case, 
the tight-binding models for the two spin components are 





12 


not necessarily consistent W e.g. assuming a mean-field 
Hubbard model, see Ref. [^, where only the on-site po¬ 
tentials can be spin-dependent. However, the obtained 
majority, minority and spinless tight-binding models for 
the ISi-lN defect are reasonable in a sense that the pa¬ 
rameters are qualitatively similar, and the values of the 
spinless parameters are somewhere between the values of 
the majority and minority parameters. Furthermore, it 
seems that at least the C-Si hopping needs to be slightly 
spin-dependent in order to obtain a qualitatively good fit 
of the band structures. 


C. Randomly positioned defects 

The tight-binding models constructed above can be 
used to study large systems containing millions of atoms, 
much larger than what can be treated by state-of-the-art 
first principles simulations. The large system sizes are 
needed to eliminate finite-size effects. Below, we study 
the electronic and transport properties of supercells with 
2000 X 1000 atoms in total and with randomly placed 
defects at a defect concentration of 0.5%. Moreover, we 
take ensemble averages over several defect configurations 
to ensure that the results are converged in the large sys¬ 
tem limit. It is, however, sufficient to average over only 
roughly five defect configuration realizations, as the sys¬ 
tems considered here are large enough for much of the 
random fluctuations to average out. 

The DOS of such large systems with randomly placed 
ISi, ISi-lN, 2Si and 2Si-lN defects are shown in Fig. [6l 
There is a strong similarity to the DOS of a periodic 
super cell with one defect presented in Figs. 0] and [5l 
However, averaging over random defect configurations 
smoothens the individual peaks of the periodic case, and 
one obtains the characteristic features of the DOS. Still, 
of course, the band structure of the periodic supercell 
and the DOS of a large system are closely related: when¬ 
ever there are relatively fiat impurity bands, there are 
impurity states at the corresponding energies, manifest¬ 
ing themselves as humps in the DOS. 

It is surprising that the DOS of a system containing 
0.5% of 1 Si defects shows only a very weak deviation from 
the DOS of pristine graphene, whereas both the majority 
and minority spin channels of the ISi-lN defect exhibit 
clear peaks. Namely, the resonant impurity states of the 
majority (minority) spin reside on the occupied (unoccu¬ 
pied) side, at negative (positive) energies. Similarly, the 
resonant impurity states of the 2Si defects correspond to 
positive energies. The 2Si-lN defect DOS in Fig. [6fb) is 
again only weakly perturbed from the pristine graphene 
DOS, but it is clearly doped by the nitrogen substitu¬ 
tions. It seems that the spiky DOS of the periodic defect 
case, as shown in Fig. [5fb), smoothens to a slight eleva¬ 
tion around E = 0 eV. 



Energy (eV) 



Energy (eV) 


FIG. 6. (Color online). Tight-binding density of states for a 
2000 X 1000 atom system with 0.5% random defect concen¬ 
tration. 

D. Electronic transport 

The tight-binding models also enabled us to perform 
transport calculations with large system sizes. We used 
the real-space Kubo-Greenwood method (RSKG)^'^— 
which scales linearly with respect to the total number 
of atoms and allows for the simulation of truly two- 
dimensional graphene with many randomly distributed 
defects. Our implementation is significantly accelerated 
by using graphics processing units (GPUs)^. 

The RSKG method can be used to efficiently compute 
the intrinsic conductivity cr(F^,t) as a function of the 
Fermi energy E and correlation time t. In most cases, 
the maximum of a{E, t) over t represents a good approx¬ 
imation of the semiclassical conductivity asc{E). The 
semiclassical conductivity is the conductivity of a sys¬ 
tem where quantum-mechanical interference phenomena 
leading to weak and strong localization are suppressed. 
The elastic mean free path le{E) can be then calculated 
using the Einstein relation for diffusive transport 

a,,{E) = ^e^p{E)v{E)UE), 


(4) 
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where v{E) is the Fermi velocity and p{E) is the total 
DOS. 

In systems with some given defect type, both the semi- 
classical conductivity and the elastic mean free path de¬ 
pend on the defect concentration n. Here, n was set to 
match the silicon concentration, and it is fixed to 0.5%. 
The defects were positioned and oriented randomly in the 
graphene lattice, as would be expected in real materials. 

The defects act as resonant scatterers^ that limit the 
charge carrier conduction around the resonant energies, 
resulting in local minima in the semiclassical conductiv¬ 
ity. The semiclassical conductivities for systems with 
many randomly positioned ISi and ISi-lN defects are 
presented in Fig. [Tfa) and with many randomly posi¬ 
tioned 2Si and 2Si-lN defects in Fig. [D^d). The sys¬ 
tem with ISi defects is again closest to pure graphene, 
showing a weak minimum only around E" = 0.8 eV. Fur¬ 
thermore, the peak at E = —0.64 eV should be even 
more singular than Fig. [3^a) shows, because in our sim¬ 
ulation the intrinsic conductivity was still increasing at 
the largest correlation times considered. For the spin- 
polarized ISi-lN defects. Fig. [3^a) shows the semiclassi¬ 
cal conductivity separately for the majority and minority 
tight-binding models. Using the same model to describe 
all defects corresponds in this case to the ferromagnetic 
ordering of the defect spin moments that would point 
along the direction of an external field. Both spin chan¬ 
nels have smaller semiclassical conductivities at energies 
that match the peaks in their respective DOS. Then, sim¬ 
ilarly to the DOS, the conductivities of the spin channels 
are completely different around zero energy, suggesting 
that ferromagnetically ordered ISi-IN defects would lead 
to distinct transport properties that depend on the spin 
states of the charge carriers. The 2Si and 2Si-IN defects 
exhibit clearly smaller semiclassical conductivities than 
the ISi and ISi-IN defects, and the minima in conduc¬ 
tivity are again found at energies that match the peaks 
in the DOS. 

From the semiclassical conductivity, we could calculate 
the elastic mean free path le{E)^ shown in Fig. [3^b) and 
(e). For the ISi case, the smallest values of le are around 
a few nanometers while the system with ISi-IN defects 
exhibits values that are even smaller than 1 nm for a 
single spin channel. Thus the ISi-IN defect is a strong 
scatterer of the charge carriers. Interestingly, the differ¬ 
ence between the spin channels is very large around zero 
energy, being around three orders of magnitude. The sys¬ 
tems with ISi defects, on the other hand, exhibit elastic 
mean free paths that exceed one micrometer, which is a 
surprising result for such a high defect concentration of 
0.5%. The elastic mean free paths in systems with 2Si 
and 2Si-IN defects are rather short in the energy win¬ 
dow shown in Fig. [H^e), as expected based on the low 
conductivities. 

To quantify the transport properties with any small 
value of defect concentration^, we calculated the scatter¬ 
ing cross section. After le{E) is obtained for a given con¬ 
centration n, the effective scattering cross section with 


dimension of [length] in two dimensions is evaluated as^ 


where rid = is the 2D number density of the de¬ 

fects, and a is the carbon-carbon distance. For relatively 
small n, which is often the case in real experiments, s{E) 
is largely independent of n. 

The scattering cross sections s{E) are shown in 
Fig. Hie) and (f). The large values of the scattering 
cross section for one spin channel in systems with ISi- 
IN defects demonstrate further that the scattering by 
this defect type is highly spin-dependent. In general, the 
scattering cross section shows that the silicon and silicon- 
nitrogen defects scatter charge carriers in clearly distinct 
ways. The ISi defects cause barely any scattering below 
E = 0 eV and scatter moderately at around E = 0.84 
eV, whereas the I Si-IN defects scatter the majority spin 
charge carriers very strongly at E = 0 eV and the minor¬ 
ity spin charge carriers moderately at E = 0.42 eV. In 
the 2Si case, scattering is very strong between E = 0 eV 
and E = 0.4 eV where, in turn, the 2Si-IN defects cause 
only little scattering. 

So far we have shown the transport results for ferro¬ 
magnetically ordered I Si-IN defects, where all the impu¬ 
rity spin moments point in the same direction. In this 
case, the charge carrier spin basis can be chosen to align 
with the impurity spins, and the majority and minor¬ 
ity spin components can be solved independently. On 
the other hand, when the external field is weak, the ISi- 
IN defects prefer antiferromagnetic ordering. However, 
the interaction energies between defect spin moments are 
rather small, roughly tens of meV, and the temperatures 
in experiments might not be low enough to manifest clear 
magnetic ordering, leading to random orientations of the 
impurity spin moments. The charge carrier spin basis 
cannot be simultaneously aligned with many random im¬ 
purity spin moments, which can lead to processes where 
the spin of the charge carrier flips in a scattering event. 
These might be relevant for spin relaxation in graphene^. 

Our RSKG implementation does not allow impurity 
spins to point in random directions, but we were able to 
simulate the case where the defects are positioned ran¬ 
domly in the graphene lattice, and each impurity spin is 
randomly taken to point either up or down. This corre¬ 
sponds to assigning the majority or the minority model 
for each defect randomly. In a sense, this can be called 
antiferromagnetic ordering of the defects, even though 
antiferromagnetism is not enforced locally. The scatter¬ 
ing cross section in this case with mixed I Si-IN majority 
and minority defects, labelled as “mixed”, is shown in 
Fig. [HI Scattering occurs now at both resonances from 
the majority and minority DOS, as one would naturally 
expect. Furthermore, the results show that one can aver¬ 
age the scattering cross sections of the two spin channels 
to obtain a total scattering cross section that is in good 
agreement with the “mixed” case. This is a special case 
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FIG. 7. (Color online). Electronic transport results for the ISi, ISi-lN with ferromagnetic ordering, 2Si and 2Si-lN defects at 
concentration of 0.5%. (a) and (d) Semiclassical conductivity, (b) and (e) Elastic mean free path, (c) and (f) Scattering cross 
section. 


of Mathiessen’s rule, with equal defect densities of two 
defect types. 

The defects considered here break the graphene sub¬ 
lattice symmetry. Therefore pseudospin is not conserved, 
and backscattering is allowed. Localization effects are 
thus expected to occur at low temperatures in large sys¬ 
tems with sufficiently large defect concentrations. The 
conductivities for the ISi-N and 2Si defects at = 0 eV 
and at = 0.15 eV, respectively, shown in Fig.[9l^a), are 
small and decrease as a function of correlation time. This 
indicates the presence of strong localization around the 
charge neutrality point (CNP), where the semiclassical 
conductivities are close to the so-called “minimum con¬ 
ductivity” cTmin = 4e^/7r/i. Converting the correlation 
time to length^, and omitting the initial region of in¬ 
creasing conductivity that corresponds to the transition 
from ballistic to diffusive transport, we obtain the length 
dependence shown in Fig.[9fb). Clearly, the conductivity 
decays exponentially as cr oc where the 2D local¬ 

ization lengths ^ can be estimated to be 92 nm for the 
ISi-lN defects and 10 nm for the 2Si defects at a defect 
concentration of 0.5%. 

In systems with ISi and 2Si-lN defects, the semiclassi¬ 
cal conductivity is large at E^ = 0 eV and at E = 0.3 eV, 
respectively, as seen in Fig. [9fa), which implies that the 
corresponding 2D localization lengths, which are expo¬ 


nentially proportional to the semiclassical conductivity^, 
are also large, presumably much larger than the phase 
coherence lengths in real systems. Therefore strong lo¬ 
calization is not expected to occur. However, for the 2Si- 
IN defects we can see weak localization, characterized by 
the slight decrease of the conductivity as a function of the 
correlation time, which is most clear around E = —0.1 
eV (not shown). At this point, the estimated 2D lo¬ 
calization length is of the order of 1000 nm at a 0.5% 
defect concentration, and it would be even larger with a 
smaller defect concentration. In conclusion, the localiza¬ 
tion effects do not play an important role in systems with 
ISi and 2Si-lN defects in real physical systems, whereas 
strong localization around the CNP is predominant for 
the 2Si and ISi-IN defect types. 


CONCLUSIONS 

By systematically optimizing the geometry of vari¬ 
ous defects in graphene containing silicon, nitrogen, oxy¬ 
gen, and hydrogen atoms, we have identified the defect 
types with the lowest formation energies. Namely, silicon 
and nitrogen prefer to fill a mo no vacancy in graphene, 
whereas oxygen and hydrogen clearly favor to stay on 
top of the graphene plane as adatoms. In general, various 
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FIG. 8. (Color online). The scattering cross section for the 
spin-polarized ISi-lN defect. The results labelled as majority 
and minority use only the corresponding tight-binding model 
at a defect concentration of 0.5%. The result labelled as mixed 
has the defects randomly chosen to be described by either the 
majority or the minority tight-binding models. The average 
is the mean of the majority and minority scattering cross 
sections, according to Mathiessen’s rule. 

impurities are effectively attracting each other to mini¬ 
mize the overall distortion of the graphene backbone. In 
this sense, nitrogen, oxygen, and hydrogen impurites and 
adatoms are trapped to silicon impurities. 

To study the electronic properties of silicon and silicon- 
nitrogen impurities, we have focused on defects where 
a silicon atom fills a monovacancy (ISi) or a divacancy 
(2Si). Energetically the most reasonable way to dope the 
ISi and 2Si defects by nitrogen is to substitute a nearest- 
neighbour carbon atom of the silicon atom with a nitro¬ 
gen atom, denoted as ISi-IN and 2Si-lN. Such nitrogen 
doping changes the electronic properties of the defects re¬ 
markably, most notably resulting in defect-induced finite 
spin moments in the case of ISi-IN defects. Each ISi-IN 
defect has a spin moment of 1.0//^ that is mainly local¬ 
ized at the neighbourhood of the silicon atom. The spin 
moments of the ISi-IN defects interact with each other 
preferring antiferromagnetic ordering. It also turns out 
that hydrogen adatoms quench any finite spin moments 
of the silicon-nitrogen defects, further enabling the use 
of such defects in spintronic devices. 

We have derived tight-binding models to describe the 
ISi, 2Si, ISi-IN, and 2Si-lN defects in graphene, and to 
evaluate electronic transport properties in realistic sys¬ 
tems with millions of atoms containing many randomly 
placed impurities. The effective scattering cross sections 
have been calculated to describe the intrinsic transport 
properties of systems with each of the defect types. Sys¬ 
tems with ISi and 2Si-IN defects have similar transport 
properties to pristine graphene, whereas the ISi-IN and 
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FIG. 9. (Golor online). Gonductivity as a function of (a) 
correlation time, and (b) length, at E = 0.0 eV for ISi and 
ISi-IN, and at E = 0.15 eV for 2Si, and at E = 0.3 eV for 
2Si-lN. The lines in (b) are exponential fits a oc e~^^^ . 


2Si defects scatter charge carriers heavily close to charge 
neutrality. Furthermore, the spin-polarized ISi-IN de¬ 
fects have strongly spin-dependent electronic transport 
properties. Namely, there are resonant impurity states 
localized at the defects in the majority spin component, 
whereas the minority spin component is a semimetal with 
similar properties to pristine graphene close to charge 
neutrality. Therefore the spin transport properties can 
be controlled by applying an external field that can tune 
the magnetic ordering of the defects. Furthermore, in 
realistically sized graphene samples, the 2Si and ISi-IN 
defects are expected to cause strong localization, whereas 
systems with the energetically preferred ISi and 2Si-IN 
defects exhibit so large localization lengths that localiza¬ 
tion effects are unlikely to be observed experimentally. 

Our results provide microscopic insights into the elec¬ 
tronic structure and transport of graphene with Si im¬ 
purities and Si-N/O/H complex defects. It is possible to 
extend graphene functionalities by tuning its properties 
through controllable introduction of defects either during 
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the growth or by using post-synthesis methods. Further 
exposure of silicon impurities to particular gases should 
give rise to other defects. On the other hand, knowing 
the signatures of particular types of defects in Si-doped 
graphene, one can use these systems in gas sensor 
or to design and further improve graphene-based nanode¬ 
vices. 
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